Protection Strategy of Sensitive Body Organs in Radiation Therapy 
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In this paper, we investigate protection strategies of sensitive body anatomy against the irradiation 
to the cancerous moving tumors in intensity modulated radiation therapy. Inspired by optimization 
techniques developed in statistical physics and dynamical systems, we deploy a method based on 
variational principles and formulate an efficient genetic algorithm which enable us to search for 
global minima in a complex landscape of irradiation dose delivered to the radiosensitive organs 
at risk. We take advantage of the internal motion of body anatomy during radiation therapy to 
reduce the unintentional delivery of the radiation to sensitive organs. We show that the accurate 
optimization of the control parameters, compare to the conventional IMRT and widely used delivery 
based on static anatomy assumption, leads to a significant reduction of the dose delivered to the 
organs at risk. 

PACS numbers: 87.56.-v,87.55.-x,87.64.Aa 



Introduction: Optimality rules in nature play crucial 
role in characterizing the dynamical pathway of complex 
systems such as enzymatic reactions [lj, protein folding 
Q , DNA repair mechanisms and carcinogenesis Q as the 
trajectories of the constituent atoms and molecules of 
these systems follow the principle of least action 0] . Op- 
timization techniques developed in recent years in statis- 
tical and computational physics have shown to be also 
useful tools in facilitating treatment delivery methods, 
in particular in intensity modulated radiation therapy 
(IMRT) d i, 0, H S, ^ IMRT is an inverse plan- 
ning approach and delivery that aims at shaping dose 
distributions better to anatomical features allowing to in- 
crease exposure of tumors and decrease the toxicity of the 
treatments @, @, 0, H, @| • Engineering treatment delivery 
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methods and the control of the leaf-velocity 
and intensity-dose-rate [1 61 ] for moving and deforming 
body anatomy [14| during hypo fractionation and IMRT 
have opened up a broad and relatively unexplored fron- 
tier for research in computational physics with applica- 
tion to radiation oncology. 

The unavoidable stochastic nature of body anatomy 
motion in living systems is one of the challenging prob- 
lems in radiation therapy and needs to be integrated as 
a part of the state of the art treatment planning sys- 
tems. Although IMRT is very important technological 
achievement, however, it faces many challanges when it 
is planned for the delivery to the body anatomies that ex- 
perience motion, particularly for treatments of thoracic 
and abdominal sites. In the presence of motion the dose 
planned and dose delivered do not match, and thus the 
clinical results can be worse than pre-IMRT treatments. 
To prevent the danger of creating secondary cancer in 
normal tissues due to this problem a proactive approach 
to utilize the organ motions has to be taken into account 
to lower the dose to sensitive structures without compro- 
mising the dose delivered to the target [HI, UH . Diverse 
criteria for quality of treatment can be utilized, some in- 



creasing the efficiency of treatment and minimizing over- 



all radiation exposure [12|, [13|, [14j while other decreasing 



cumulative dose distribution over org ans at risk (OAR) 



to minimize treatment toxicity [15l . |161 | . The recently pro- 
posed algorithms aim to change the paradigm of IMRT 
to moving body anatomy. Taking into account the fact 
that multiple delivery strategies are possible to moving 
body anatomy, the goal is to use the additional degrees 
of freedom in the dynamical system of body anatomy 
and multi-leaf collimator (MLC) systems, resulting from 
motions of body organs, to the advantage of treatment. 
Thus our goal is not to achieve dose distributions in mov- 
ing body anatomy as good as it is possible for static 
anatomy but be able to deliver dose distributions to mov- 
ing body that are advantageous to any dose distribution 
possible for static anatomy. The solution strategies can 
lead to complex optimization processes, many support- 
ing diverse landscapes of multiple local minima. Among 
few possible procedures leading to achieving global solu- 
tions, or satisfactory local solutions, we explore explicitly 
dynamical pathways that breed offspring control param- 
eters that efficiently converge to advantageous solutions. 
These so called genetic operations [T3], such as repro- 
duction, suppression, mutation, and the mating crossover 
lead to efficient solutions of our multi-parameter dynamic 
treatment optimization problems. 

We present a strategy based on controlling the velocity 
of multi-leafs and intensity dose rate of beams in dynamic 
multi-leaf collimator (DMLC) deliveries that allow to re- 
distribute beneficially the irradiation to the target and 
to organs at risk over subsequent phases of the body mo- 
tion. To simplify our computations, we assume in this 
study that motion of the body anatomy is known be- 
fore delivery, though this assumption is not compulsory 



for our approach [15l . |16| | . In particular by using the ve- 
locity of the leaves and the intensity dose rate as the 
control parameters in optimization process we propose a 
strategy that the dose delivered to the organ at risk can 
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be reduced significantly relative to the delivery based on 
static body anatomy assumption. 

Method: In the following we present a method based 
on optimization algorithm that mixes minimization of the 
time of delivery (and thus minimization of the contami- 
nating irradiation) and minimization of dose delivered to 
OAR for exposure of moving body anatomy. We use ve- 
locity of the leaves and the intensity beam-dose rate (or 
intensity rate) a as control parameters in optimization. 
We consider the motion of the target to be aligned with 
the movements of MLC leaves in one-dimension. 

Optimizing time of delivery: 



Refs. [12 

vestigated the mathematical model for the time-optimal 
control of a single leaf pair (or a set of multiple inde- 
pendent leaf pairs) for the case of moving rigid target. 
They explored deliveries for time-optimal DMLC strate- 
gies subjected to constant intensity rate. To preserve the 
integrity of the delivered intensity in these solutions, the 
following leaf velocity (vf), the leading leaf velocity (vl), 
and the target velocity (vt) have to satisfy the relation 
given below 



vf(x, mp) 



vl(x, m L ) - VT(m L ) 
1 + r(x)[v L (x,m L ) - v T (m L )} 



+ v T (m F ).(l) 



In Eq. |T]), mi and mp denote the time measured in cu- 
mulative monitor unit (MU) at which the leading and fol- 
lowing leaves pass over a point on the target denoted by x 
(defined in the target frame of reference), I'{x) = dl/dx, 
and mp — niL = I{x). Here I{x) is the planned inten- 
sity, intendend to be delivered to the target. The veloc- 
ities are expressed in cm/MU, i.e., v\ = dx\/dm where 
A £ {T, L, F} and x l if) is the position of the edge of 
the leading (following) leaf in the laboratory frame of 
reference, xt is the position of the target. Eq. ([1]) is in- 
variant under variation of a. Because of the mechanical 
design of MLC, leaves can not move faster than a maxi- 
mum accessible velocity v max . This condition enforces a 
lower-bound m m i n for the final moment of delivery m/ 
such that rrif > m m j n . rn m i n decreases with increas- 
ing Umax- To achieve a unique solution corresponding 
to m/ = m m - m with constant a, one of the leaves moves 
with the maximum velocity v max , while the other moves 
with appropriately matched velocity 12j . It follows 



dxF 
drriF 



dx L 
dmr, 



'"max— «t(™l) 



for case 2 
for case 1 



(2) 



Umax" VT ("If) 

l-i'(*)[« m »-«T(m F )] + V T( m L) 



for case 1 
for case 2 



(3) 



Here v max is expressed in MU. Trajectories of the leaves 
can then be easily calculated by integrating the veloci- 
ties. Conditions defined as case-1 and case-2 denote the 
appropriate subsets of the Cartesian product of spatial 



variables spanning the whole range of target in the mo- 
tion, and the interval of time (exceeding the time of de- 
livery). We refer to Refs. [ij, EE 13 for details. Ap- 
plying the transformation dm/dt — a or equivalently 
m = J Q dt'a(t r ) in Ecl([T]) leads to leaf velocities expressed 
in physical time t [1 61 ] 



VL(x,t L ) 



v F {x,t F ) _ 

"l + I'(x)[^ 



VT(t L ) 



(**) 



VT{tp) 

a{tp) 



(4) 



where tup) is the physical moment at which leading (fol- 
lowing) leaf passes over point x on the target. 

Optimizing irradiation dose delivered to the organ at 
risk: Ref. [l5| investigated an algorithm based on a 4D 
DMLC leaf sequencing for the purpose of minimizing 
the irradiation dose delivered to OAR. The optimization 
has been carried out by segmentation of the velocities of 
leaves in case when delivery was characterized bycon- 
stant a. According to the results presented in Ref.[l_5|, it 
is possible to reduce the delivered dose to OAR to 80% 
of the dose minimized relative to the best starting point 
of treatment delivery during all phases of motion of the 
treatment anatomy. Here we adopt similar strategy as 
presented in Ref. [151]. However with the inclusion of a 
as a new controllable parameter, we illustrate that the 
reduction of the dose delivered to OAR can be improved 
significantly. 

We formulate the problem of minimizing the irradia- 
tion dose to the OAR by the variational principle [4[ . We 
define the integral of irradiation dose delivered to OAR as 
an objective function (or cost-function which is required 
to be minimized) 



Doar[xl,xf,xt,u] = / 



dx Iqar(x) 



dx (mp — mi,){x"). (5) 



Here x denotes a point on OAR (in laboratory frame of 
reference), and £Ioar € [xi,X2] defines the effective area 
and boundaries of OAR. Formally Doar is a functional 
of xl, xf, xt and a. Among these parameters we can 
vary xl and a [as the motion of the target is not exter- 
nally controllable, and xp is assumed to be a function 
of xl, given by Eq. [TJ or Eqs. (J2j{3j)] . Geometrically 
Eq. ([5]) represents the area surrounded by the leaf's tra- 
jectories and the boundaries of OAR, equivalent to the 
area of OAR under radiation. The optimal solution in 
form of system of Eulcr-Lagrange equations [J] can be 
determined via 



5D 



OAR 



= 0. 



(6) 



The functional variation has been carried out with re- 
spect to independent variables xl, and a. Minimizing 
Dqar requires longer time of delivery with respect to the 
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treatment plans in which minimization of dose delivered 
to OAR is not taken into account. Because analytical 
solutions of this problem are not achievable, a compu- 
tational method based on genetic algorithm is employed 

S3' 

Computational aspects: The numerical optimization is 
carried out in two steps. In first step we start with an 
initial constant a, and calculate the optimal leaf trajec- 
tories defined by velocities vl, and Vp, using Eqs.©, 
and ([3]). For the optimization of velocities we divide the 
target into -ZV scg segments. In each segment we assign 
a random number to and calculate Vp. We use ge- 
netic algorithm that iteratively transforms a population 
of (vl,vf) with fixed-length N seg each with an associ- 
ated fitness value defined as Dq\ r , into a new popula- 
tion of offspring (vl,vf) using the Darwinian principle 
of natural selection and using operations that are pat- 
terned after naturally occurring genetic operations, such 
as mating crossover, suppression and mutation [171 ] . As 
a result, optimal velocities corresponding to constant a 
are identified. 

In second step we use leaf trajectories obtained in the 
first step, together with the constant beam dose rate 
as assumed initially, as a starting strategy denoted by 
(a, vl,vf). In this step we calculate the optimal vari- 
able dose rate with appropriate leaf trajectories. The 
result of this calculation lead to an updated strategy 
plan (cx,Vl>Vf). The connection between strategies in 
steps one and two can be formulated by an appropriate 
transformation which keeps Eq. [4] invariant, and pre- 
serves the integrity of the dose delivered to target. In 
the target frame of reference this transformation is given 
by mt = rriL, and rhp — rap. Using MU to physical time 
transformation, we find an equivalent transformation for 
the velocities of the leaves 

\ [v L (t L ) ~ »r(t£)] = - [v L {t L ) - v T (t L )} ■ (7) 
a a 

Here v L {th) = dx L {t L )jdt h , and v L (t L ) = dx L (t L )/dt L - 
From Eq.([7]) one can deduce 

Xh{ih) - x T (h) = x L (t L ) - x T (t L ). (8) 

Similar equations holds for Vp, and xp. Similar to the 
first step, at point x on the target, delivered dose to 
the target is calculated by I(x) — rhp — rh^. Here 

™>L{F) = Jq L<f) dt'a(t') where the leading (following) leaf 
arrives at x on the target at t L ^ F y Note that the trajecto- 
ries calculated in second step are required to preserve the 
integrity of I(x). These equations imply that the trajec- 
tories of the leaves in target frame of reference, expressed 
in MU, are invariant under variation of the intensity dose 
rate. This in turn means that shape of I(x) (that is deter- 
mined by the mutual relationship between leaf velocities) 
is kept invariant (integrity of delivery is preserved). The 
procedure described above, base on steps 1 and 2, can be 
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FIG. 1: The double parabola intensity used in the calcu- 
lations. The solid line indicates the intensity map intended 
to be delivered and the line with stars shows the intensity 
actually delivered to the moving target. 



repeated in a loop that continuous till it converges to a 
satisfactory solution. 

Numerical results: For illustration of the method we 
consider a moving rigid target with length xj = 4 and 
a static OAR in laboratory frame of reference. For sim- 
plicity, the trajectory of the target is assumed to be one- 
dimensional and is given by xt = sin 4 (0.6i + 1.8). The 
dimension of OAR is assumed to be 2.5 < x < 3.5. The 
length and time are expressed in cm and second respec- 
tively. The planned intensity to be delivered to the tar- 
get is a double parabola as illustrated by solid line in Fig. 
[TJ We also consider the maximum accessible velocity of 
the leaves and intensity dose rate as f ma x = 3 cm/sec, 
and a max = 1000 MU/min. For the optimization of the 
leaf velocities, we divide the target up to seven segments 
(iV se g = 1 to 7). In each segment we associate ten ran- 
dom velocities to the leading leaf subjected to vl < t> max . 
The number of possible combinations for and thus the 
size of search space is given by lO^ 8 and exceeds up to 
10 7 . For the optimization of a, a step- wise beam dose 
rate with one step has been assumed with four parame- 
ters: the initial and final moments, and the intensity dose 
rates of the pulse. The optimal value of these parame- 
ters are calculated by the random search optimization 
procedure. The optimized delivery plans shows that the 
reduction of the dose delivered to OAR is possible to be 
as low as 45% for N scg = 3, 5, and 6. For illustration 
we present the numerical results in Figs. [5] and [3J cal- 
culated based on the optimization of vl with N scg = 5 
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FIG. 2: Trajectories of the leaves shown in the laboratory 
frame of reference expressed as a functions of cumulative num- 
ber of monitor units m. The delivery strategy is based on 
optimized time of delivery and irradiation dose delivered to 
the organ at risk. Figure (a) shows trajectories of leading 
and following leaves calculated based on uniform intensity 
dose rate a — 10 MU/sec (dashed line), and a — 2 MU/sec 
for 3.5 < t < 4.5 sec and a = 10 MU/sec elsewhere (solid 
line). The latter has been calculated based on optimization 
of step- wise dose rate that minimizes intensity to static organ 
at risk. Both trajectories are calculated based on unoptimizcd 
leaf velocities. We observe that modulation of the intensity 
dose rate lowers the irradiation dose delivered to ORA. In fig- 
ure (b) optimization of the leaf velocities has been performed 
by segmentation of the target into five pieces. In each piece 
the optimal vl, and vf have been calculated using genetic al- 
gorithm, with a — 10 (dashed line), and optimized intensity 
dose rate (solid line). In the latter, the self-consistent opti- 
mum intensity dose rate is a = 13 MU/sec for < t < 2.5 
sec and a — 10 MU/sec for t > 2.5 sec. A comparison with a 
strategy based on uniform intensity dose rate a — 10 MU/sec 
with unoptimized leaf velocities shown in figure (a) by dashed 
line, and the self-consistent optimized leaf velocities and in- 
tensity dose rate shown in figure (b) by solid line, shows that 
the dose delivered to the organ at risk (the area surrounded 
by the leaves trajectories and the boundaries of the OAR) can 
be reduced significantly. 



FIG. 3: Profiles of intensity delivered to static organ at risk 
calculated for static target (green dashed line), optimized with 
respect to beam dose rate (black and red dashed lines), and 
optimized with respect to both of leaf velocities and beam 
dose rate delivery (purple solid line). The integral of dose 
delivered to the organ at risk is calculated for each case and 
is given by 12.8, 9.5, 8, and 6 MU. We observe that the opti- 
mization based on leaf velocities and intensity dose rate may 
reduce the irradiation dose delivered to organ at risk down to 
45%. 



corresponding to optimal a given by 

f 13 MU/sec < t < 2.5 , , 

a ~ \ 10 MU/sec t > 2.5 ^ ' 

Conclusion: In this work we showed that the internal 
motion of body anatomy during therapy can be used as 
a valuable information in protecting the sensitive organs. 
To minimize the error in IMRT delivery to moving targets 
we developed strategies based on controlling the velocity 
of multi-leafs and intensity rate of beams that allow the 
irradiation of the target to be properly re-distributed over 
subsequent phases of target motion. This general ap- 
proach utilizes multiple solutions that provide the same 
dose to target while optimizing these delivery parameters 
that can improve the quality of the therapy. Our calcu- 
lation based on optimization of the control parameters 
using genetic algorithm shows that the dose delivered to 
the organ at risk can be reduced significantly compare to 
the delivery based on static anatomy assumption. 
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